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We introduce a multiscale Monte Carlo algorithm to simulate dense simple fluids. The probability 
of an update follows a power law distribution in its length scale. The collective motion of clusters 
of particles requires generalization of the Metropolis update rule to impose detailed balance. We 
apply the method to the simulation of a Lennard- Jones fluid and show improvements in efficiency 
over conventional Monte Carlo and molecular dynamics, eliminating hydrodynamic slowing down. 



Both molecular dynamics and Monte Carlo are widely 
used in the study of fluids. The first, which is often pre- 
ferred for large scale studies, allows one to access both 
static and dynamic information; however use of a large 
time step leads to systematic errors in the thermodynam- 
ics Q . Monte Carlo has the advantage of converging to 
the exact equilibrium distribution without any system- 
atic errors; recently the introduction of collective updates 
has led to substantial improvements in the efficiency of 
Monte Carlo simulation of fluids. Particularly notable is 
a cluster algorithm for hard spheres Q and its generaliza- 
tion Q to Lennard-Jones fluids. These papers introduced 
updates based on reflection or rotation of groups of par- 
ticles with respect to a randomly chosen centre. The al- 
gorithms are particularly well adapted to the simulation 
of dilute, heterogeneous systems. They work less well 
for dense homogeneous fluids. In this Letter we intro- 
duce collective updates which simulate such fluids more 
efficiently. 

Our method displaces blocks of particles on a scale, 
I, which is intermediate between the dimensions of the 
particles and the simulation cell. By moving groups of 
particles we accelerate the dynamics of long wavelength 
fluctuations. Three technical problems are to be solved 
in the implementation of the method: Firstly we must 
move particles without substantially increasing the en- 
ergy of the system; otherwise the Metropolis algorithm 
will refuse the trial. Secondly we create a set of moves 
that are reversible, in order to apply detailed balance. 
Thirdly we calculate the Jacobian that is implied by our 
collective moves, which complicates the counting of the 
number of states. We treat each of these three points 
before presenting an implementation and studying the 
relaxation of long wavelength density fluctuations. 

Rigidly displacing a block of particles a distance e gen- 
erates a mismatch in the structure of the fluid at the in- 
terface between the moved and stationary particles. This 
mismatch rapidly leads to a high energy cost and low effi- 
ciency. We mitigate the problem by introducing updates 
in which the mismatch is spread over the whole volume 
i d in d dimensions. We choose updates along a single 
coordinate direction a 
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for particles, i within i of cq. The continuous function g 
has the properties g(0) = 1, and g(r) = for |r| > 1. e 
is a random amplitude, cq a randomly chosen centre in 
the simulation volume. Thus only particles within I of 
the point Co move. 

The proposed displacement, eg generates a map x — > 
G(x) on the configuration space of the system, where x 
designates the N x d coordinates of the N particles. The 
simplest manner of generating dynamics that converge 
to the Boltzmann distribution is the imposition of de- 
tailed balance between pairs of states that are linked in 
both directions by a transformation. However as yet G is 
not paired with an inverse transformation, G _1 . We con- 
struct explicitly the inverse. We add to the moves eq. Q 
a second set found by solving eq. 0) for rt as a function 
of r[. We then choose between direct and inverted moves 
with probability 0.5. 

We now turn to the imposition of detail balance. Stan- 
dard derivations consider a discrete space: In the pres- 
ence of transitions between two states i, j with energies 
Et and Ej one requires that 

p(i -> j)7ti =p(j -> i)7tj (2) 

where the occupation probabilities are given by the Boltz- 
mann weight 7tt = e — /Z. One solution for the transi- 
tion rates p(i — > j) is the Metropolis choice: 
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In the continuum one should count the number states 
in some neighborhood dx of x and not just work with 
the probability density at x. Consider two states x, x' 



linked by the transformations x' = G( 



x' = G 
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The transformation distorts the the volumes dx, and dx' 

at 



which are related by the Jacobian Jg = TjTT ' P a * r 
{G,G~'} bring two volumes dx and dx/Je into corre- 
spondence. The number of states near x is 7t x dx; near x' 
there are 7t X 'dx/jG. The transition probabilities should 
then satisfy 

p(x->x')7tx=p(x'->x)^ (4) 

Jg 

in order to generate the probability density n. This equa- 
tion has a solution 
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with T = (3 _1 the temperature. We note that Jacobian 
weighting factors are used in Monte Carlo simulations 
of polymers Q when working with torsional degrees of 
freedom. Standard continuum Monte Carlo updates of 



the form r{ 



e have Jg = 1 and fall within our 



formalism. 

The hnal choices that need to be made concern e and I. 
The energy cost of a trial is estimated from the strain in- 
duced in the fluid by the transformation: S ~ e Vg(r/£) ~ 

We 
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e/i. This energy varies as J S 2 ~ S 
match this to the thermal energy and choose 
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setting the scale of the random amplitude e |l2( . We sam- 
ple I from a distribution: Standard arguments give the 
density of states in Fourier space as N(q)dq ~ q d_1 dq. 
When this is expressed as a function of a length scale 
I = 1/q we find N{i)di = d£/£ d+1 . We cut off the dis- 
tribution at short distances I < l c and sample I with 
probability 



P(£) 



l c < I < L/2 



Fluctuations at all length scales are then sampled at the 
same rate. 

The asymptotic cost of the algorithm per trial can be 
estimated by recognizing that the effort needed to up- 
date the scale I is 0(l d ) when interactions are short 
ranged. The average work per collective update is then 
$l d P{l)&l - l d log(L/l c ), where L is the system size. The 
effort needed to perform the collective updates diverges 
only weakly with system size. 

We implemented the algorithm, simulating Lennard- 
Jones particles with the potential 



U = 4e 
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truncating the interactions at a radius of 2.5cr and shift- 
ing U by a constant to eliminate the discontinuity 
We use length units of cr and energy units of £lj- Tem- 
perature is measured in units of ei_j. Periodic boundary 
conditions were imposed. 
A convenient choice for g is 



g[ r /t) = 1 
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for r < t, g = otherwise. This allows the direct in- 
version of eq. QJ. However, we also implemented more 
elaborate expressions for which analytic inversion was not 
possible. For these cases we performed a numerical in- 
version of eq. Q using Newton iteration. Even here gen- 
erating the trial moves was a minor contribution to the 
total computational cost of the algorithm. 



The Jacobian is calculated as a product of individual 
particle contributions. For a displacement in the direc- 
tion a 
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The first simulations were performed on a very small sys- 
tem in order to check the correctness of our arguments 
as to the need for a Jacobian in the Metropolis criterion. 
To do this we multiplied the term in log Jg by a numer- 
ical prefactor y. y = 1 corresponds to the algorithm as 
described above, y — corresponds to neglecting the dif- 
ference between the volumes dx and dx'. We took five 
particles in a two dimensional box of side L = 10 and 
simulated using a conventional Monte Carlo algorithm 
as well as the multiscale algorithm. As can be seen from 
Figure n use of the correct Jacobian is essential. These 
trials on small systems allowed the generation of high 
statistics runs, from which we checked the correctness 
of the code using both the energy and partial structure 
factors. 
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FIG. 1: Variation of the mean energy, (U), as a function of the 
prefactor of logic- The correct average energy is generated 
with y = 1 . Errors result if one neglects the Jacobian, y = 
0. Solid line: local Monte Carlo. Points: multiscale Monte 
Carlo. Statistical errors small than symbol size. Vertical line 
guide to eye for y — 1 . 

In order to study the efficiency of the algorithm we 
simulated a two-dimensional fluid with cell size L = 60, 
80, 120 and 180 at a density p = 0.755, T = 0.5. We 
performed simulations with conventional, single particle 
Monte Carlo, molecular dynamics and multiscale Monte 
Carlo. In the simulations we measured the static struc- 
ture factor S q = | Y.i e lqri | . This correlation func- 
tion is used to calculate the compressibility for q — > 0. 
In charged systems similar, long wavelength, polarization 
correlation functions must be measured with high accu- 
racy in order to deduce the dielectric constant @|. S q 
is slow to equilibrate at small q due to the locality of 
particle motion in conventional Monte Carlo: A hydro- 
dynamic description of the dynamics involves a local con- 
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served quantity, the density, which diffuses: i.e. Model 
B dynamics |t|. Relaxation rates vary as T q = Dq 2 , 
with D the collective diffusion coefficient. Similarly in 
molecular dynamics propagating sound waves imply that 

r q 2 = cV- 

Slow long wavelength fluctuations have consequences 
for the equilibration of the energy. For Monte Carlo the 
autocorrelation of the energy for the mode q can be ex- 
pressed in the form C q (t) = C(Dq 2 t)withC(0) =0(T 2 ), 
so that in d dimensions the two time energy-energy cor- 
relation function is give n by C u = J C(Dq 2 t) d d q. One 
finds C u (t) ~ 1/t d/2 |fj|]. It is the integral of Cu. that de- 
termines the convergence rate of the average energy || . 
In two dimensions, where J C u (t) dt diverges for large t, 
short ranged functions converge slowly due to hydrody- 
namics. Good thermodynamic statistics require that all 
modes in the sample are equilibrated. The situation is 
better in three dimensions, the integral is dominated by 
short wavelength fluctuations. 

We verified that the algorithm worked correctly when 
only using updates from the multiscale distribution. 
However this turns out to be sub-optimal for the perfor- 
mance of the simulation. Standard Monte Carlo methods 
work very well at the scale of single particles; it is only 
at larger distances that they become slow. We stochasti- 
cally mixed a standard Monte Carlo algorithm with the 
multiscale method using a lower cut off l c = 7.5. We 
used a proportion of 1000 local updates for each multi- 
scale update. With these choices of parameters the exe- 
cution time of the program increases approximately 30% 
for L = 40 and 60% for L = 1 80 compared with the purely 
local code. Acceptance rates of all moves were adjusted 
to be approximately 40%. We verified that the accep- 
tance rate of multiscale moves was almost independent 
of I by binning acceptance. We determined the integrated 
autocorrelation times of S q with a blocking algorithm || 
and plot the inverse time, T q in Figure (0). As expected 
local Monte Carlo leads to slow relaxation of long wave- 
length modes, with P q = Dq 2 . 

For the molecular dynamics simulations we used a 
time step X = 0.006 and a damping coefficient for the 
Langevin thermostat of 0.2. The time step is typical 
for low accuracy molecular dynamics studies, but larger 
than that which must be used for accurate studies of 
thermodynamic properties, for example T = .0014 in 
The dynamics (in Fourier space) of a system of particles 
coupled to an external thermostat lead to a dispersion 
law familiar from the damped harmonic oscillator [l4j : 
Kq 2 — tncu— pcu 2 = withr| a friction coefficient, K a com- 
pressibility and p the mass density. For large q modes 
are under-damped so that w 2 — Kq 2 /p, corresponding to 
propagating sound waves. At small q we crossover to an 
over-damped "Brownian" regime where iv =iKq 2 /r|. In 
our simulations we placed the crossover near q 2 = 10~ 2 
by performing several trial runs with different r\. In this 
way we ensure propagative behavior over most of the 



range of the curve. At short length scales molecular dy- 
namics is slower than local Monte Carlo, however it does 
a better job of relaxing long wavelength correlations. We 
adjusted the number of time steps in the molecular dy- 
namics simulation so that the total simulation time was 
identical to the conventional Monte Carlo code. 

The multiscale results are again displayed with a time 
scale of inverse sweeps, for a "clock time" comparison one 
must lower the values of V by a factor between 1 .3 and 
1 .6 depending on L. The multiscale algorithm displays 
two regimes. At short wavelengths local Monte Carlo is 
active leading to standard q 2 behavior. At long wave- 
lengths the algorithm eliminates hydrodynamic slowing 
down- all long wavelength modes relax at very similar 
rates. For L = 1 80 the multiscale algorithm relaxes long 
wavelength modes 40 times faster (clock time) than con- 
ventional local Monte Carlo. Multiscale is also faster 
than molecular dynamics at this scale. 
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FIG. 2: Inverse integrated autocorrelation time in sweeps, 
T q , for S q . Conventional Monte Carlo: +. Multiscale: □. 
Molecular dynamics: A. Conventional algorithm has slow, 
diffusive modes. Multiscale has decay rate independent of q 
at long wavelengths. Dashed line: r ~ q 2 ; dot-dashed line 
r ~ q. Time scale for Monte Carlo inverse sweeps. System 
sizes from L = 40 to L = 180. Simulations of 2 20 « 10 6 
sweeps, taking up to two weeks of simulation time for the 
largest systems. 

Given the success of the algorithm we tried to find 
other updates which couple even more strongly to density 
fluctuations. We implemented a trial update in which 
motion is purely radial about a random centre Co. If r 
denotes the position of a particle with respect to cq wc 
tried 

r' = r+er(1-T/«) (II) 
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for r/i < 1 , r' = r otherwise. The results were disap- 
pointing. Similar mediocre results were found for purely 
rotational updates of the form 



e(1 -r 2 /£ 2 



(12) 



with the angular position of a particle seen from Co. 
Mixtures of eq. (|TT)l and eq. (fT^fl were just as poor. 
Clearly, not all multiscale updates are created equal. 

Whilst we performed the most detailed simulations in 
two dimensions we did implement the three dimensional 
version of the code and used it to simulate a single system 
with L = 40. We performed both conventional and mul- 
tiscale Monte Carlo and found, Figure qualitatively 
similar results to Figure |2 It is to be noted that three 
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FIG. 3: Inverse integrated autocorrelation time in sweeps, 
T q , for S q in three dimensions. Conventional Monte Carlo: 
+ . Multiscale: □. T = 1 , L = 40, l c = 2.5. Line F ~ q 2 . 
Fq tends to a constant at long wavelengths for the multiscale 
algorithm. 

dimensional simulations require considerably more com- 
puting resources for a given system size. The number of 
particles is larger, as well as the number of neighboring 
particles within 2.5cr. 



In this Letter we have introduced a wider choice of 
Monte Carlo moves for the simulation of simple fluids. 
A good choice of updates leads to increases in the effi- 
ciency of simulations in a manner reminiscent of multi- 
grid Monte Carlo simulation of lattice models . Using 
conventional Monte Carlo we demonstrated that the re- 
laxation rate varies as T = Dq 2 ; density fluctuations in 
large systems take 0(L 2 ) sweeps to equilibrate. Mul- 
tiscale updates, with an average cost of V~ log L/l c per 
update, eliminate long wavelength slowing down, leading 
to an algorithm which is asymptotically faster. 

Our method permits more general implementations of 
preferential sampling in which updates are performed 
more often in the neighborhood of an interesting site or 
surface. It may also have application in heterogeneous 
colloidal systems where co-motion of small and large par- 
ticles or molecules may be an alternative to cluster meth- 
ods at high densities. 
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On repeating the argument with a block of rigidly dis- 
placed particles we find the exponent d — 2 replaced by 
d — 1 . Large blocks can not move as far. 
We also checked this decay in two dimensions using 
Kawasaki dynamics for a lattice gas. 
We neglect internal friction of the fluid which leads to 
weaker dampening. 



